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Planetary-scale giant storms erupt on Saturn quasiperiodically. There have been at least six recorded occurrenc- 
es of past eruptions, and the most recent one was in 2010, with its whole life span captured by the Cassini 
mission. In 2015, we used the Very Large Array to probe the deep response of Saturn’s troposphere to the 
giant storms. In addition to the remnant effect of the storm in 2010, we have found long-lasting signatures 
of all mid-latitude giant storms, a mixture of equatorial storms up to hundreds of years old, and potentially 
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an unreported older storm at 70°N. We derive an ammonia anomaly map that shows an extended meridional 
migration of the storm’s aftermath and vertical transport of ammonia vapor by storm dynamics. Intriguingly, the 
last storm in 2010 splits into two distinct components that propagate in opposite meridional directions, leaving 


a gap at 43°N planetographic latitude. 


INTRODUCTION 

Every 20 to 30 years, a giant storm erupts on Saturn, creating enor- 
mous cloud disturbances with a clear head marching forward until 
it wraps around the whole planet (1, 2). The most recent eruption 
was in 2010, with visible cloud activities lasting for more than 6 
months (3-5). At that time, the Cassini spacecraft was still orbiting 
Saturn and the on-board instruments, including the visible camera, 
radio instruments, infrared camera [composite infrared spectrome- 
ter (CIRS)], and Cassini Radar, returned exquisite details of the as- 
tronomical feast. The next event is expected to happen within 10 to 
20 years, based on prior statistics and the understanding of the 
storm dynamics (1, 2, 6). In the stratosphere, more than 10 K tem- 
perature anomalies appeared in a concentrated stratospheric beacon 
and lasted for 2 to 3 years after the storm (7, 8), causing a substantial 
change in the chemical composition in a local area (9, 10). 

Despite similarities between the giant planets, storms of these 
sizes are unique to Saturn. One peculiar aspect of the last giant 
storm on Saturn is its “dehydration” effect, referring to the observa- 
tion that the storm depletes the condensable vapors in the atmo- 
sphere, in particular, the ammonia vapor. Since ammonia vapor 
provides sufficient microwave opacity at the altitudes probed by 
the Cassini/Radar instrument, the observed high brightness tem- 
peratures in the Cassini Radar’s radiometric observation at 2 cm 
was interpreted as being caused by a completely dry atmosphere, 
i.e., no ammonia vapor, down to at least ~2 bar (11, 12). 

This observation was puzzling because moist convection would 
bring ammonia-rich air from deep levels to shallower levels, which 
should enrich the ammonia vapor concentration in the shallower 
layers rather than deplete it. The apparent contradiction can be ex- 
plained by the interstorm meridional circulation. Unlike dry con- 
vection, moist convection is asymmetric with respect to updrafts 
and downdrafts: The initial updrafts carry moist and saturated air 
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upward, but the subsequent downdrafts push dry and unsaturated 
air downward. The process may oscillate several times, leaving 
behind an overall unsaturated atmosphere (6). Similar arguments 
have been used to explain a depletion of ammonia on Jupiter in 
the upper versus lower atmospheric layers (13). 

Another possible mechanism to remove ammonia vapor is by 
forming ammonia-rich mushballs (14, 15). Protected by a thin ice 
shell, ammonia-rich (up to 30%) liquid (16) can be sealed within 
centimeter-sized hail-like precipitation and be effectively transport- 
ed to the deep atmosphere. Sequestration of ammonia vapor into a 
cold ammonia-water solution at shallower layers and release of it at 
deep layers where the temperature is high may explain the depletion 
of ammonia gas. This mechanism also indicates colder but 
ammonia-rich deep layers due to evaporation. A warm and dry 
upper layer sitting on top of a cold and moist deep layer stabilizes 
the atmosphere and prevents the occurrence of moist convection for 
many years. However, the post-storm stable stratification would be 
eroded away by the infrared cooling of the atmosphere and by tur- 
bulent mixing that homogenizes the ammonia concentration. Both 
processes are slow because the cooling flux is low and the atmo- 
sphere is nearly inviscid. 

We examine the long-term evolution of the atmosphere of 
Saturn after the 2010 storm using the upgraded Karl G. Jansky 
Very Large Array (VLA). No big storms were observed after the 
last giant eruption, and the post-storm evolution of the atmosphere 
would be governed by the aforementioned slow process: infrared 
cooling and mixing. Saturn's atmosphere is a huge heat reservoir, 
even for the shallow part at altitudes above the 30-bar pressure 
level. Using the emitted power, 5 W/m”, 5 years of cooling will 
only lead to a change in the column-averaged temperature by 
about 0.2 K (17). Thus, we suspect that observed variations in 
brightness temperature, while sensitive to both temperature and 
ammonia variations, are predominantly due to variations in the 
ammonia vapor. 


RESULTS 

VLA observations 

VLA observations of Saturn were obtained on two dates: 24 to 25 
January 2015 when the array was in its compact configuration 
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(C&B), and 28 to 29 May 2015 with the array in its more extended 
configuration (B&A). These are hybrid configurations in which the 
antennas on the east and west arms have already been moved to the 
more compact configuration, but those on the north arm remain 
extended. Such array configurations are ideal to observe sources 
at relatively low declinations. In January (C&B), the planet was ob- 
served in the Q (39.8 to 47.6 GHz), K (18 to 26 GHz), U (12 to 17 
GHz), and X bands (8 to 11.8 GHz), and in May (B&A) in S (2 to 2.5 
GHz), C (4 to 4.2 GHz), X (8 to 11.2 GHz), and U (12 to 17 
GHz) bands. 

These data were presented previously in the analysis of Saturn's 
ring system (18), with a summary of observations provided in their 
table 1. All data were initially processed and calibrated using the in- 
ternal VLA calibration pipeline, followed by using the MIRIAD 
software to self-calibrate, clean, and map the data (19). On the 
basis of a preliminary model fitting of Saturn's disk (18), we deter- 
mined the disk-averaged brightness temperatures in each band, as 
summarized in table S1. 

Figure 1 shows current and previous observations of Saturn's 
disk-averaged radio emission, together with a model spectrum for 
comparison. This disk-averaged model assumes a dry adiabatic at- 
mosphere and a uniform distribution of ammonia vapor up to the 
cloud base, referred to as the homogeneous isentropic model. Above 
the cloud base, the ammonia vapor follows its saturation value 
uniquely determined by the temperature. We use the Voyager 
radio occultation profile (20) for temperatures at altitudes above 
the 1-bar pressure level and use an adiabatic extrapolation for tem- 
peratures below the 1-bar level. The VLA channels do not probe 
deep enough to see through the condensation level of water, so 
we used a nominal value of 10 times the solar water abundance 
(21), estimated from the dynamics of Saturn's giant storm (6). 
The disk-averaged brightness temperatures of the C and S bands 


are sensitive to the ammonia concentration above the water cloud 
level near 20 bars. The x’ fitting of both bands reduces to X? < 1 for 
an ammonia abundance of around 200 parts per million by volume 
(ppmv), which corresponds to an enrichment factor of 1.5 com- 
pared to the solar N abundance (see the Supplementary Materials, 
section A, for a study on the ammonia concentration). However, 
this value may not be truly representative of Saturn's deep 
ammonia abundance because the S band is not long enough to dis- 
tinguish various deep ammonia abundances [see the comparison of 
two other models using 150 and 400 ppm deep ammonia abun- 
dance (22, 23)]. Nonetheless, Saturn's disk-averaged emission can 
be fitted adequately by a homogeneous isentropic model, while Ju- 
piter’s atmosphere shows substantial depletion of ammonia vapor 
in the upper layers compared to the deep atmosphere (24-26). 

Figure 2 shows images of Saturn in different projections for all 
six radio bands with the ring contribution removed. The left column 
displays residual disk images as viewed on the sky with a limb-dark- 
ened disk subtracted (based on disk-averaged parameters as shown 
in table S1); the middle column projects the residual disk image into 
a cylindrical projection, with rotational smearing of about 60° in 
longitude; the right column zooms in on the polar region north- 
ward of 50°N. We obtain the limb of the disk by tracing the sharp 
edge of the full disk image and overlaying that on the residual disk 
image. We focus our analysis on the northern hemisphere where the 
line of sight was free of ring materials. 

Clear signatures of bright and dark latitudinal bands are found 
from the residual disk images at wavelengths shorter than 10 cm (C, 
X, U, K, and Q bands). At the S band, the banded structure is barely 
visible except for the shadow of Saturn’s ring in the southern hemi- 
sphere. However, this may be due to the low spatial resolution of the 
S-band image (~10° in latitude). The C-band image starts to show 
one bright latitudinal band between 20°N and 40°N and a relatively 
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Fig. 1. Disk-averaged observation and model of Saturn. (Left) Saturn's disk-averaged brightness temperature as a function of wavelength. VLA observations from this 
paper are shown by the orange points. The blue dots are older data (22, 38, 42-44). Wilkinson microwave anisotropy probe (WMAP) (45) and Planck (46, 47) data are shown 
in chocolate and magenta colors, respectively. All data have been corrected for cosmic microwave background radiation, and use the full Planck equation (rather than the 
Rayleigh-Jeans approximation) to convert flux densities to brightness temperatures. (Middle) Contribution function of the radio bands and the dry adiabatic temperature 
profile. (Right) Ammonia and water profiles. 
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Fig. 2. VLA maps of Saturn in radio bands. (Left column) Disk-removed images of Saturn at S, C, X, U, K, and Q bands sequentially from top to bottom. The beam pattern 
for each band is drawn in the lower left corner of the disk image as ellipses. The orange circle in each panel indicates the limb of the disk. (Middle column) Cylindrical 
projection of the residual brightness temperature. (Right column) Polar projection of the residual brightness temperature. The lowest latitude is at 50°N in the projection. 
These images are rotationally smeared; the center column shows the approximate longitudinal range covered by each observation. A longitude-resolved map for the 
combined U and X band data was shown in (48). B/W color scale shows the residual brightness temperature with a limb-darkened disk subtracted. 
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bright equator. The dark band in the southern hemisphere is due to 
the ring occulting the emission of the planet. 

The X-band image features an extremely bright latitude band 
between 40°N and 50°N. The brightness temperature of this latitude 
band is about 15 K warmer than the neighboring bands, showing 
that ammonia vapor must be substantially depleted in the atmo- 
sphere. Looking at the polar projection map in the right column, 
we also find a bright ring surrounding the north pole, which we at- 
tribute to the deep signature of Saturn's hexagon, although the hex- 
agonal shape is unclear at the resolution of this image. 

The highest spatial resolution is achieved in the U-band images, 
in which the bright 40°N and 50°N band breaks into two distant 
narrow bands with one at 40°N and the other at 44°N. This dual- 
band structure is only present in the U-band images and is realistic 
as it consistently shows up at all longitudes in the cylindrical pro- 
jection map. The bright ring-like structure surrounding the pole 
near 78°N now takes a hexagonal shape in the U-band image. 
Since its position is colocated with the visible image of the 
hexagon, we interpret it as being the deep root of Saturn's polar 
hexagon. At a higher altitude (~100 mbar), the Cassini/CIRS in- 
strument found warmer temperatures (~4 K) than the neighboring 
latitudes (27). This pattern resembles our finding at the U band. Al- 
though the altitude probed by Cassini/CIRS is higher than in our 
radio map, it is likely that the hexagonal jet remains warm in tem- 
perature and extends deep (at least ~2 bars) into the lower 
troposphere. 

The 44°N bright feature in the U-band image also exists in the K- 
band image, but the 40°N bright feature disappeared, probably due 
to a decrease in spatial resolution. The bright equatorial zone is the 
most prominent feature in the K-band image, accompanied by a 
broad dark band centered near 20°N. The Q-band image is noisy 
and the banded structure is almost indistinguishable, although at 
these frequencies we probe similar pressure levels to that of the 
U band. 


Zonal-mean residual emission 

Figure 3 shows the zonal-mean latitudinal scans of the radio emis- 
sion across spectral bands, after a limb-darkened disk was subtract- 
ed from each image (Fig. 2). The noises are illustrated as gray ticks in 
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each band. Since the emission angle changes substantially across the 
disk, we average over 60° in longitude around the central meridian 
of the disk image. For the S band, we increase the width of the lat- 
itude bin to obtain at least eight pixels to have a robust estimate of 
the noise. Thus, the uncertainty reflects the noise in the radio map, 
which is much smaller than the calibration uncertainty. Any longi- 
tudinally dependent features are smeared out because these scans 
are zonally averaged. But since these scans cover several spectral 
bands, they reveal changes in brightness temperatures over a 
range of altitudes. 

These zonal-mean scans at the different wavelengths form, in 
essence, a spectral data cube; we evaluate these at the emission 
angle between the local zenith and an Earth's observer. Multiple 
bright features, such as the 40°N and 44°N bands, are discernible 
from the spectral scans and can be correlated across spectral space 
from top to bottom. In general, a feature appears bright (high 
brightness temperature) when the ammonia vapor, the main micro- 
wave absorber, is depleted compared to the average, or the temper- 
ature is high. For example, the brightness temperatures are higher at 
45°N in the Q, K, U, and X bands but lower than average in the C 
band compared to the neighboring latitudes, suggesting two end- 
member scenarios: (i) The ammonia concentration is depleted at 
higher altitudes (sensed by Q, K, U, and X bands) and enriched at 
lower altitudes (sensed by C band), or (ii) the temperature is higher 
at higher altitudes and lower at lower altitudes compared to the 
neighboring latitudes. The VLA observation alone cannot resolve 
the ambiguity, but these two scenarios are always related: Evapora- 
tion/sublimation of ammonia rain droplets/ice crystals would si- 
multaneously cool the atmospheric temperature and enrich the 
ammonia vapor concentration. 

To have a quantitative assessment of the physical properties of 
Saturn's atmosphere, we perform a beam deconvolution followed 
by a differential spectral inversion to recover the ammonia concen- 
tration anomaly (see Materials and Methods). Fitting the brightness 
temperature anomaly rather than its absolute value bypasses the 3 to 
5% (28) calibration uncertainty and leverages the precision of the 
measurements. Figure 4 shows the deconvolved spectral scans in 
orange lines and the width of the Gaussian beam in red bars. De- 
convolution has the largest effects on K and U band data due to the 
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Fig. 3. Residual zonal-mean brightness temperatures at six bands. The mean is obtained by averaging over 60° in longitude around the central meridian of the 
observed disk. Gray ticks in each panel reflect the 10 uncertainty estimates from all longitudes within each latitude bin. The spacing of the step segment represents 


the pixel size of the disk image. 
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Fig. 4. Deconvolved brightness temperature anomaly. From (top) to (bottom), the radio bands are ordered as Q, K, U, X, C, and S. Orange lines show the deconvolved 
brightness temperature anomaly. Cross-hatched patterns indicate the noise level in a latitude bin. Green lines are reconvolved brightness temperature anomalies using 
the red lines for each band. Red bars at the lower left corner show the FWHM of the Gaussian beam. Anomaly is with respect to the mean emission of the northern 
hemisphere. 
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high precision, small beam width, and large spatial variation of the 
brightness temperatures. After the deconvolution, the peak ampli- 
tude of the brightness temperature anomaly is raised and the phase 
(position of the brightness temperature peaks) of the spectrum is 
shifted. Deconvolution has a relatively small effect on Q, X, C, 
and S bands given their relatively smooth variation and/or large 
beam size. 

The variation of the brightness temperatures in Fig. 4 reflects the 
variability of Saturn's atmospheric composition or temperature on 
top of an imprecisely determined baseline. From Fig. 1, we find that 
the brightness temperatures of a homogeneous isentropic model at- 
mosphere could adequately explain the disk-averaged brightness 
temperatures of Saturn. Therefore, we choose the baseline model 
of the differential fitting to be the simple homogeneous isentropic 
model. This choice also provides an additional benefit such that the 
anomaly map can be compared both vertically (altitude) and hori- 
zontally (latitude) since the baseline model has a uniform distribu- 
tion of ammonia and is isentropic. We further fix the temperature 
profile and only examine the distribution of ammonia because the 
observed brightness temperature anomaly is far greater than the in- 
crease in temperature due to the giant storm in 2010 (2.4 K) (29) or 
the latitudinal variability of temperature at the 0.5-bar pressure level 
(~2 K) measured by the Cassini/CIRS instrument (30) based on the 
collision-induced continuum absorption of hydrogen at far infra- 
red. The Cassini/Visible and Infrared Mapping Spectrometer in- 
strument provides somewhat deeper observation of thermal 
emission at near infrared, but deriving the meridional variability 
of temperature is impossible because of the degeneracy with PH; 
and aerosols (31). The formula and details of the differential 
fitting are provided in Materials and Methods. 
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Anomalous ammonia concentration map 

Figure 5 shows the ammonia anomaly map with Saturn's zonal 
wind profile in System HI (32) for comparison. We find three dis- 
tinct patches of ammonia anomalies, between 0°N and 12°N, 
between 25°N and 45°N, and between 55°N and 72°N. There are 
six occurrences of giant storms on Saturn in the past 150 years: 
the 1876, 1933, and 1990 storms erupted between 0°N and 12°N; 
the 1903 and 2010 storms erupted between 35°N and 42°N; and 
the 1960 storm erupted near 58°N. We find ammonia anomaly sig- 
natures for all mid-latitude eruptions and an equatorial mixture 
with their long-time evolution. Specifically, we find that the 
ammonia anomaly propagates meridionally in the increasing direc- 
tion of the zonal wind, meaning that the meridional velocity is as- 
sociated with the vorticity of the flow. The propagation of the 
ammonia anomaly is usually stopped at the position of another 
zonal jet. On the basis of these observations, we suspect that there 
could be another storm that has occurred near 70°N latitude that 
was before 1876, the first recorded storm. 

Another finding is a division between atmospheric layers at ~3 
to 5 bars caused by the giant storms: At altitudes above 3 to 5 bars, 
the ammonia vapor is depleted, and below that level, the ammonia 
abundance is enriched. During the eruption of the 2010 storm, the 
Cassini/Radar observed highly depleted ammonia gas at roughly the 
1- to 2-bar level (12). But where did the ammonia vapor go? If 
ammonia is hidden within the ammonia-rich precipitation (14), it 
would reevaporate and enrich the ammonia concentration in the 
deep atmosphere. This is exactly what we have seen using longer- 
wavelength radio observations: After a giant storm, ammonia gas 
is indeed transported from the upper atmosphere to deeper layers. 
The layering of the atmosphere indicates stable stratification, and 
the simplest form is that the atmosphere is divided into two distinct 
layers. On the basis of the observation of the layering, we estimate 
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Fig. 5. Ammonia vapor concentration anomaly with respect to a homogeneous model. The zonal wind profile is indicated by the green curve at the top, with its scale 
along the right y axis. The locations and the years of the previous giant storms are annotated in the figure. Storms in 1876, 1933, 1903, and 1960 are denoted by their 
center locations. The equatorial storms (1876, 1933, and 1990) (49) and the 2010 storm (3-5) are displayed as a latitude range. 
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the convective mass flux of the 2010 storm using an energy balance 
calculation. 

The last 2010 giant storm increased Saturn's emitted power in 
the stormy band by ~9.2%, which translates into a 2.4 K increase 
in physical temperature (29). From the ammonia anomaly map, 
we find that the transition pressure level is at about pp ~ 3 bars. 
Therefore, we posit that at least 3 bars of atmospheric mass is 
warmed up either by the direct thermodynamic response to the 
latent heat release or by secondary circulation such as the compen- 
sating subsidence through adiabatic heating (6). The amount of 
latent heat released per unit area per unit time is equivalent to the 
product of the mass flux of the upward motion and the mass mixing 
ratio of water, multiplied by the latent heat of water. Equating that to 
the increase of column temperature yields 


car ees = Lm,q,St (1) 
£ 

where the heat capacity of the hydrogen-helium atmosphere is c, = 
1.3 x 10*J/ (kg K); the increase in temperature is AT = 2.4 K. The 
upward convective mass flux is denoted as m,,. The spatial area 
covered by the storm is S = 3 x 10'° m? based on the 300,000 km 
circumference of the planet and the 10,000 km width in latitude, 
and t = 6 months is the duration of the storm. Saturn's gravity at 
40°N is g = 11 m/s*. The latent heat of water is L = 2.5 x 10° J/kg, 
and the mass mixing ratio of water participating in the convection is 
qw. If we use the deep water mixing ratio, which is about 10% of the 
total atmospheric mass (6), the convective mass flux is 


E cpATPp 
 qugtL 

Assuming that the detrainment level of the upward convective 
mass flux occupies one scale height at the physical temperature of 
100 K, H = RT/g = 34 km, where the atmospheric pressure is p = 0.4 


bar and the atmosphere density is p = p/(RT) = 0.1 kg/m°. The as- 
sociated cloud expansion rate of the giant storm in 2010 is thus 


= 2.2 x 10-* kg/(m’s) 


(2) 


u 


. ins 
A = = 1.9 x 10? km?/s 
pH 


This value is very close to the measured cloud expansion rate of 
A = 200 km?/s based on the sequence of cloud images over the 
whole lifetime of the storm (3, 4). The division levels of other 
storms are located at lower altitudes than that of the 2010 storm, 
indicating potential subsidence over the years following a giant 
storm. The subsidence could be driven by radiative cooling and 
convective adjustment. For example, the downward propagation 
of a stable interface separating a convective interior and a stable 
upper atmosphere was simulated by a one-dimensional convective 
adjustment model (6). In addition, as described in the theory, the 
downward propagation is stopped at the water condensation level, 
which is around 20 bars assuming a 10 times solar water abundance. 

The third finding revealed by the ammonia anomaly map is that 
a mixing barrier seems to form where the zonal wind velocity 
crosses zero. This is most evident near the two zero-crossings of 
the zonal wind profile flanking the 40°N westward jet. Specifically, 
the post-eruption evolution of the 2010 storm is intriguingly split 
into two distinct components. The northern one propagates north- 
ward and the southern one propagates southward, leaving a gap 
in between. 


(3) 
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The gap is located at 43 + 0.4°N based on the high-resolution U- 
band image of Saturn. The location is very close to the zero-velocity 
crossing at 41.8°N using the System III rotation period (10 hours, 39 
min, 24 s) determined by Voyager's measurement of Saturn's kilo- 
metric radiation in 1978 (32). Saturn's rotational period is uncertain 
because Saturn's dipole magnetic field is not tilted relative to its ro- 
tation axis (33). Other indirect measurements based on atmospheric 
dynamics, response to gravity, and ring seismology generally place 
Saturn's rotation period in between 10 hours 32 min 35 s and 10 
hours 34 min 13 s (34-36), which corresponds to a zero-crossing 
latitude near 44.2°N. It is unknown what causes the separation 
and the formation of the gap. We raise a conjecture that the gap 
could indicate where the zonal wind is zero. This can be tested by 
a numerical model studying the long-term evolution of the post- 
eruption tropospheric dynamics. 


DISCUSSION 

For most of the time, Saturn's atmosphere looks hazy and feature- 
less to the naked eye in contrast to Jupiter's colorful and vibrant at- 
mosphere. This picture changes when we look at Saturn using a 
radio eye. Despite bland looking at the visible wavelengths, different 
latitudes on Saturn show dramatic contrast in the radio emission. In 
2015, we found substantial brightness temperature variations— 
greater than the temperature increase caused by the 2010 storm— 
as a function of latitude and radio wavelength based on spatially re- 
solved radio images of Saturn, as obtained from the VLA. Via a 
spectral inversion study, we created a latitude-altitude map of 
ammonia concentration anomalies from these images and correlate 
the anomalies to the occurrences of past giant storms in the north- 
ern hemisphere. The equatorial mixture of storms and storms at 
mid-latitudes are featured by a depletion of ammonia at altitudes 
above the 3- to 5-bar level and an enrichment of ammonia below 
that pressure level. The most likely scenario is that ammonia is 
transported from the upper atmosphere to the lower atmosphere 
via precipitation and reevaporation. The ammonia concentration 
anomaly also seems to propagate meridionally across latitudes 
during long-term post-eruption evolution. Specifically, the past 
2010 storm is divided into two distinct components, with one prop- 
agating southward and the other northward, leaving a gap at 43°N. 
These findings call for more investigations into the large-scale at- 
mospheric dynamics during the post-eruption state of Saturn's at- 
mosphere. Our result shows that Saturn's tropospheric dynamics 
may be substantially different from Jupiter's. The Juno microwave 
radiometer revealed a correlation of brightness temperature 
anomaly with zones and belts. However, on Saturn, the brightness 
temperature anomalies at radio wavelengths are dominated by 
giant storms. 

The observation we took in 2015 was unable to view Saturn's 
southern hemisphere due to its rings. We are looking forward to 
another observing opportunity in 2025 when Saturn's rings will 
be viewed edge-on and both hemispheres can be observed from 
Earth. Since there are no detections of giant storms in the southern 
hemisphere, we expect that brightness temperature variations are 
smaller in the southern hemisphere than in the north. 


7 of 12 


EZOZ ‘81 ISNBNY UO BIO‘aOUSTOS: MMAM//:Sd}IY woy popeoqumoq 


SCIENCE ADVANCES | RESEARCH ARTICLE 


MATERIALS AND METHODS 

VLA data acquisition and calibration 

The VLA data were acquired in January and May 2015 and present- 
ed previously with the analysis of Saturn's ring system, referred to 
here as Zhang et al. (18). As discussed in Zhang et al. (18), all data 
were initially processed and calibrated using the internal VLA cal- 
ibration pipeline. They subsequently used the MIRIAD software to 
self-calibrate, clean, and map the data (19). Zhang et al. (18) com- 
bined all eight channels in each spectral window, as well as several 
spectral windows to increase the signal-to-noise ratio (SNR). In U 
and X bands, the final bandwidth per window was 1 GHz; in the K 
band, it was 4 GHz; in the Q band, all windows were combined with 
a resulting bandwidth of almost 8 GHz. In C and S bands, all spec- 
tral windows were combined as well, resulting in a bandwidth of 
0.25 GHz in each band. For X and U bands, data were obtained 
in two different array configurations. To increase the sensitivity 
for atmospheric study, we combined the data obtained on different 
dates by scaling the flux density of the source to match both obser- 
vations. This flux density is integrated over the solid angle extended 
by the radio source, and thus captures the changes in the change of 
Saturn's diameters on different dates. 

The flux density scale was calibrated on the radio source 3C286, 
with the standard VLA flux calibrator scale (28, 37). The radio 
source J1558-1409 was used as the complex gain calibrator 
(phases). Internal and absolute uncertainties in the flux densities 
are believed to be better than ~3% at most frequencies and 5% at 
K and Q bands. We adopt these, somewhat conservative, uncertain- 
ties in our work during the fit to the disk-averaged brightness tem- 
peratures. These uncertainties only affect the precision of the 
baseline inversion and do not affect the accuracy of the retrieved 
ammonia anomalies through the differential spectral inversion. 

Zhang et al. (18) used a “planet + ring" model to self-calibrate the 
data and determine the best-fit parameters for both the planet and 
Saturn's rings. The planet's brightness temperature T, of the unoc- 
culted portion of the planet's disk was modeled as T, = To + Ty(Up); 
where To and T; are model parameters of the uniform and limb- 
darkened disk, respectively, and u, is the cosine of the emission 
angle. The rings were modeled using the Simrings code (38). Mod- 
eling the rings this way enabled Zhang et al. (18) to determine the 
ring properties of Saturn's main rings, as, e.g., their optical depth, 
particle size distribution, composition, and porosity. In our work, 
we use final calibrated UV data from Zhang et al. (18) and subtract 
the model for the rings to construct the best possible data for 
Saturn's atmosphere. Rather than using the combined C&B and 
B&A array configurations for the X and U data, we instead focus 
on only the B&A configuration for these bands, although we use 
the more compact array to calibrate the flux densities. 

On the basis of Zhang et al.'s model (18) fits of Saturn's disk, we 
determined the disk-averaged brightness temperatures in each spec- 
tral window, as summarized in table S1. Since Saturn blocks out the 
cosmic background radiation (CMB), and the flux density is mea- 
sured relative to the background, we add the temperature corre- 
sponding to the CMB to the disk-averaged brightness 
temperature. We further used a disk-averaged value of 0.667 for 
the cosine of the emission angle of the disk when combining the 
To and Tı temperatures into a disk-averaged value. While we 
derive here (table S1) the disk-averaged brightness temperature in 
each spectral window, to increase the SNR in the maps, we 
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combined all spectral windows per frequency band, resulting in a 
total of six maps. 


Deconvolution algorithm 

All radio maps presented here have different resolutions, i.e., they 
are all "smoothed" by their respective beams, which differ from 
wavelength to wavelength. To obtain a quantitative assessment of 
the physical properties of Saturn's atmosphere, we need to convolve 
models of Saturn's atmosphere with their respective beams in an it- 
erative process to determine the cause of observed variations in 
brightness temperature with latitude. To speed up this process 
and to mitigate the effect of the beam picking up cold sky signals 
when close to the limb, we instead deconvolve the latitudinal 
scans in each wavelength band and then use a differential spectral 
inversion to obtain an ammonia anomaly map, i.e., the difference in 
ammonia profiles relative to the nominal (200 ppm) value as a func- 
tion of latitude that best matches the data. 

We denote the beam-convolved zonal-mean brightness temper- 
ature anomaly as T(x), where x is the projected distance normal to 
the observer's line of sight, referred to as the sky coordinate. The 
corresponding point-wise brightness temperatures are denoted as 
T,(y). We approximate the VLA beam pattern as a Gaussian func- 
tion, denoted as G(y; s), where s represents the full width at half 
maximum (FWHM) of the beam and is related to the SD of the 
Gaussian function by s = 2v 2ln20. Using the sky coordinate x, 
the beam convolution gives 


oe 100 

Fils) = | TH)G—ys)ay (4) 
To deconvolve the observed brightness temperature, we use the 

Fourier transform and regularized optimization. The Fourier trans- 

form of Eq. 4 is 


F{Ts}(w) = F{T}(w)- F{G}(w) (5) 


where w is the angular frequency and the product is performed 
element-wise at each frequency. 

In principle, one can just divide the Fourier transform of the 
convolved beam, F {Ty}, by the Fourier transform of the Gaussian 
beam, F{G}(w,), to get the Fourier transform of the pointwise 
brightness temperature, F{T,}(w;), and then invert it for the point- 
wise brightness temperature. However, this method suffers from 
numerical instabilities for high-resolution maps. On the one 
hand, high-resolution observations contain more information, but 
on the other hand, it introduces more Fourier modes to invert. 
Usually, the increase of the information content is not enough to 
constrain the additional Fourier modes, and therefore, a direct di- 
vision yields an unstable solution. We use regularized optimization 
to overcome the difficulty. In discretized frequency domain, w;, Eq. 
5 becomes 


F{T,}(wi) = F{Ty} (wi) F{G}(w;),i = 0...n (6) 


Now, we cast the equation above into a matrix multiplication 
form such that 
F {Ty} (wi) = Gir F{To} (w) (7) 


where Gj = F{G}(w;)ô; are the coefficients of a diagonal matrix. 
Instead of solving Eq. 7 exactly, we try to obtain a regularized 
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least-square solution that demands vanishing amplitude for high- 
frequency modes such that 


Gi F{To} (w) = F{T} (mi) (8) 


subject to small F{T,}(w;) for large w;. The way to do it is to add a 
penalty term for the amplitude of F{T,}(w;) as a function of its w;. 
By doing experiments with triangle functions (see the Supplemen- 
tary Materials, section B), we find that a penalty proportional to fre- 
quency squared works well for the deconvolution. Mathematically, 
the cost function of the minimization is 


n-1 
Cost function = x? + VS“ [wiF {T} w) 


i=0 


= S165 F(T) — FE Mwy? 


i=0 


(9) 


+ eS wr (To Hn)? 


where À is a tunable parameter governing the strength of the regu- 
larization. Written in matrix notation, the least-square solution to 
Eq. 9 is 


F{T,} = (GTG + VW'W) 'GTF{T,} (10) 
where W is a diagonal matrix with its elements being W; = w;ð;. 
The component form of Eq. 10 is 
F{G}(wi) z 
FIL Wj) = FIT Wi 11 
{Ty} (wi) FIG w) + ew {Ty} (wi) (11) 


Compared to Eq. 6, the A? term in the denominator of Eq. 11 
provides the stability of the solution. For low-frequency modes, 
where F{G}(w) > Aw”, Eq. 11 reduces to Eq. 6. For high-frequency 
modes, their amplitude would decay as w™°. The cutoff frequency, 
Woa occurs at F{G}(w.) = Awe. It is useful to define a cutoff wave- 
length, Le = = to denote the spatial resolution that can be recov- 
ered from the deconvolution. Spatial variations at wavelength 
smaller than L, will be suppressed by the deconvolution, and vari- 
ations greater than L, will be preserved. Note that Eq. 11 bears re- 
semblance to the Wiener deconvolution, which has been used 
extensively in the imaging processing community to sharpen a 
noisy image (39), such as restoring the Hubble Space Telescope 
(HST) images (40). However, the primary goal of the Wiener de- 
convolution is to clean up a distorted image due to random noise, 
and our goal is to recover the underlying pointwise brightness tem- 
perature from a (beam-convolved) observation. For a Gaussian- 
shaped beam, its Fourier transform is 


+00 ii x2 
F{G}(w;0) = L Lep( eww iwx)dx 


1 
= exp (- sow) 


We calculate the discrete Fourier transform of the convolved 
brightness temperature, F{T,}, using fast Fourier transform 
(FFT) at discrete frequencies w;. Then, we use Eq. 12 to calculate 
the Fourier transform of a Gaussian beam at the same discrete fre- 
quencies and apply the result to Eq. 11 to get the spectral amplitudes 


(12) 
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of the pointwise brightness temperature. Last, inverse FFT is per- 
formed to recover the deconvolved brightness temperatures. 

To examine the performance of the deconvolution, we recon- 
volve the beam using Eq. 5 and compare that to the original 
(beam-convolved) brightness temperatures. The results are dis- 
played in fig. S2. The choice of the penalty strength A? depends 
on the desired cutoff wavelength L,, which is related to the width 
of the radio beam. From the synthetic study of a triangle signal func- 
tion (Supplementary Materials, section B), we conclude that the 
spatial sampling resolution can be between one-quarter and one- 
half of the beam's FWHM. For S, C, and X bands, we use one- 
half of the beam’s FWHM as the cutoff wavelength, and for U, K, 
and Q bands, we use one-quarter of the beam’s FWHM as the cutoff 
wavelength. The effect of cutoff wavelength can be illustrated by 
looking at the distribution of the Fourier amplitudes. For 
example, fig. S2 shows the real part and the imaginary part of the 
Fourier amplitudes of the beam-convolved brightness temperatures 
(blue circles and blue crosses). The cutoff wavelength of the S-band 
data is visually identified at L. = 0.1 L, where L is the diameter of 
Saturn's disk, and that for the X band decreases to about 0.02 L. 
The beam’s FWHM at the S band is about 1/5 L, and therefore, 
any variations of wavelength smaller than 1/10 L are suppressed. 
Deconvolution does not have an effect over wavelengths greater 
than 1/10 L. Large-scale (compared to the beam width) features 
are preserved in the deconvolution. 


Differential spectral inversion method 

After the deconvolution, we perform a differential fitting using the 
algorithm below. Let Ty(y, u; A) be the brightness temperature at 
latitude y, at cosine emission angle u = cos®, and at wavelength À. 
We decompose T, into a baseline component T? subject to a cali- 
bration uncertainty and a fluctuation component T} subject to a 
measurement precision such that 


Toly; u, À) = Tp (uo, A) + (A) + TO; HA) + 6(y,A) (13) 


where [lo is the averaged emission angle, e(X) is the random variable 
representing the calibration uncertainty as a function of wavelength 
only, and 8(y, A) is the random variable representing the measure- 
ment precision, which is a function of both wavelength and latitude. 
We assume that T? >> T} and the calibration uncertainty is much 
larger than the precision of the measurement: e > ô. To represent 
the fluctuation, we request the average of T} to be zero. Therefore, 
T? + £ can be obtained by spatial average and T} + 6 is the residual 
after subtracting the mean. 
Suppose that the radiative transfer equation is written as 


Ta; tA) = fla); WAI 
where q represents the vertical profile of ammonia vapor, which is a 
function of latitude. Suppose that there exists an unknown profile, q, 
satisfying 


(14) 


F (45 A) = Th (to, A) + (A) (15) 
we use the mean value theorem to obtain 
Toly; m, A) =F(G u, X) +f [ga(y)ss Nia) -— 41 a6) 


where q,(y) is between q(y) and ĝ. Upon matching the terms in Eqs. 
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13, 15, and 16, we find that 


Fli os lay) — ål = TiO: wA) + 6(y,) (17) 


On the other hand, if Eq. 15 can be approximately fitted by a 
known profile, qo, within the calibration uncertainty such that 


Ff (do; u, A) © TH(to, A) + €(A) (18) 


Then, for any profile q(y), using the mean value theorem again, 
we have 


Fa); u, A) =F (do; u, X) + lao); u, Alla) — 40] (19) 


where q, is between qo and q(y). Note that q, may not be known, but 
f (qv) can be precisely evaluated from the radiative transfer calcula- 
tion. Subtracting Eq. 19 from Eq. 17 yields 


F'{qa(y); u Ala) — ål- Flao); u, Mla”) — qo] = 


T05 wA) — (Fla): mA] — flgo; w Al} + 80,2) 

Since qo and q are sufficiently close (up to the calibration uncer- 

tainty) and if the brightness temperature function f(q) is approxi- 

mately linear for a range of q that yields small variation of T} 
compared to T}, we can safely assume that 


F'O): A] = flg); pA] 
without knowing q and q, exactly. Denote the ammonia concentra- 
tion anomalies as 


(20) 


(21) 


5"q = q0) — qo (22) 
5°q = 40) -å (23) 
the differential fitting formula reads as 
flas): u, A] (8q — 8" q) = THO pA) — Fla): uA] 
= flgo; u, Al} + 8,4) (24) 


and the minimization function is 


Minimizeg|T4(y; u, A) — {fla(y); mA] — flgo; p, AH) 80,4)? 
(25) 


If Eq. 21 is strictly valid, then the fitted ammonia anomaly 5™q con- 
verges to the true ammonia anomaly 5“q. Otherwise, the inference 
would be biased. The validity of Eq. 21 depends on how close qo is to 
qand how linear the brightness temperature function is with respect 
to the ammonia concentration. Figure $7 shows the optical depths 
across various VLA radio bands. It shows that different bands are 
sensitive to different pressure levels. At about 10 bars, S band's 
optical depth is close to 1. At about 20 bars, S band's optical 
depth increases to 4. Linearity is best observed when the optical 
depth is close to 1. We terminate our inversion algorithm at 25 
bars, as beyond this point, the sensitivity is lost and the linearity as- 
sumption is no longer valid. Nevertheless, using the differential 
fitting formula, Eq. 25 bypasses the large calibration uncertainty € 
and leverages the measurement precision 6. 

The fitting and radiative transfer were performed using the 
Markov chain Monte Carlo (MCMC) method as described in (24, 
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41). First, we divide the atmosphere into five layers with pressure 
boundaries at 25, 12.2, 5.98, 2.93, 1.43, and 0.7 bars, respectively. 
These pressure levels are chosen on the basis of log spacing 
between 25 and 0.7 bars, which denote the lower and upper boun- 
dary of the spectral inversion domain. Second, we fix the ammonia 
concentration at 200 ppm at 25 bars (see the Supplementary Mate- 
rials, section A, for the choice) and vary the ammonia concentration 
at the other pressure levels. The ammonia profiles in between these 
pressure levels are smoothly interpolated. In total, we have six sam- 
pling variables. Third, we calculate the brightness temperatures of 
this model atmosphere at the emission angle specified by the ob- 
serving geometry for each frequency band. Last, we minimize Eq. 
25 based on the MCMC algorithm. More discussions on the inver- 
sion parameters are provided in the Supplementary Materials, 
section C. 


Supplementary Materials 
This PDF file includes: 

Sections A to C 

Figs. S1 to S8 

Table S1 
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